In this tutorial, we will see how we can perform an uncertainty analysis and a sensitivity analysis for a crop model
Based on the result of the sensitivity analysis, we will then see how we can optimize the model, to reduce the gap between the previsions and the observations
This tutorial follows this lesson, where we learned how to code a crop model with R
Introduction
As weather input, we will use the same data set as in the previous tutorial:
library(tidyverse)day_sowing<-92# Sowing after 1st Mayday_harvest<-287# Harvest ~ mid-octoberdata<-readr::read_delim("https://raw.githubusercontent.com/BjnNowak/CropModel/main/Weather_DesMoines.csv",delim=";")%>%arrange(LST_DATE)%>%mutate(day_number =row_number())%>%select(day_number,tas = T_DAILY_MEAN,rsds = SOLARAD_DAILY)%>%filter(day_number>=day_sowing, day_number<=day_harvest)# Check first lineshead(data)
Previously, we only coded the first part of the model, to determine the number of leaves for corn. In order to save some time, I wrote a script with the rest of the equations, that can be loaded with the source() function as shown below:
# Load script from githubsource('https://raw.githubusercontent.com/BjnNowak/Lessons/main/scripts/crop_model_function.R')# Take a look at the fun_model() functionfun_model
The arguments of the fun_model() function are the following:
name # Scenario name
data # Climatic variables inputs
GDD_1leaf # Phyllochron
C # Crops’s capacity to intercept light
RUE # Radiation use efficiency (gDM.MJ-1)
nthresh # Max number of leaves
Yield prediction
We will now predict the maximum potential corn yield for our climate conditions, unsing the fun_model() function:
Daily grain yield previsions (in tons) are stored in the yield_t column
The results may be plotted as follows:
p1<-ggplot(data=baseline,aes(x=day_number,y=grain_t))+geom_line()+labs( title ="Corn yield estimation for DesMoines",x ="Day number",y ="Potential max yield (t.ha-1)" )+theme_light()p1
Yield prediction
If you are only interested in the final yield, it may be extracted as follows:
# Extract final yieldbaseline%>%tail(1)%>%pull(grain_t)
[1] 13.63428
Uncertainty analysis
An uncertainty analysis is used to assess the variability of the previsions of the model
An uncertainty analysis provides a confidence interval for model predictions
To do this, the variability of all relevant parameters is considered simultaneously
Uncertainty analysis
For this case study, we will only focus on the variability of two parameters here:
C, which reflects the crops’s capacity to intercept light. It depends on the light extinction coefficient (reflecting light penetration in crop’s foliage), the individual leaf area, and the plant density. We will will consider here that C may vary between [0.06; 0.18] (ie up to 18% of light intercepted) ;
the Radiation Use Efficieny (RUE), in gDM.MJ-1, which estimates the conversion of the intercepted radiation into biomass. We will will consider here that RUE may vary between [1; 3].
Model structure
For this case study, we will only focus on the variability of two parameters here:
C, which reflects the crops’s capacity to intercept light;
the Radiation Use Efficieny (RUE), in gDM.MJ-1, which estimates the conversion of the intercepted radiation into biomass.
Uncertainty analysis
Now that we have set the variability limits for these parameters, let’s create a data table with random combinations of the parameters in their variability range
N=100# Number of drawsanalysis<-tibble(# Column with C parameter valuesC_draw=runif( N, # Set number of draws0.06, # Min value0.18# Max value ),# Column with RUE parameter valuesRUE_draw=runif( N,1,3 ),# Empty column for corresponding final yield values# (to be filled later)yield=0)
Uncertainty analysis
We fill now complete the final yield values corresponding to each {C,RUE} combination with the fun_model() function, using a for loop
for (i in1:N){ analysis$yield[i]<-fun_model(name="DesMoines", data=data, # C and RUE combination based on row indexC=analysis$C_draw[i],RUE=analysis$RUE_draw[i],# For this example, GDD and number of leaves are stableGDD_1leaf =50,nthresh =16 )%>%# Extract final yield onlytail(1)%>%pull(grain_t)}
Uncertainty analysis
We may now plot the result of our uncertainty analysis
Considering the uncertainty of the parameters C and RUE, the final yield lies within the following limits:
p2<-ggplot(analysis,aes(y=yield,x=1))+geom_boxplot()+geom_jitter()+labs( title ="Uncertainty analysis",subtitle="Each point corresponds to one simulation",x="",y ="Potential max yield (t.ha-1)" )+theme_light()p2
Sensitivity analysis
A sensitivity analysis is used to determine the effect of uncertainty about the variables or parameters of a model on the final result
A sensitivity analysis identifies the parameters that have the greatest influence on the final model result
To do this, the variability of each relevant parameters is considered separately
Sensitivity analysis
We will start by initiating data table for both parameters C and NUE
N=100# Number of drawsanalysis_C<-tibble(C_draw=runif(N, 0.06, 0.18),yield=0)analysis_RUE<-tibble(RUE_draw=runif(N,1,3),yield=0)
Sensitivity analysis
Then we will fill these tables with a for loop
for (i in1:100){# Analysis for C analysis_C$yield[i]<-fun_model(name="DesMoines", data=data, # Variable parametersC=analysis_C$C_draw[i],# Stable parametersGDD_1leaf =50,RUE=2,nthresh =16 )%>%tail(1)%>%pull(grain_t)# Analysis for RUE analysis_RUE$yield[i]<-fun_model(name="DesMoines", data=data, # Variable parametersRUE=analysis_RUE$RUE_draw[i],# Stable parametersGDD_1leaf =50,nthresh =16, C=0.12 )%>%tail(1)%>%pull(grain_t)}
Sensitivity analysis
We may now plot the result of our sensitivity analysis (using scale() to plot the results for both parameters on a relative scale)
p3<-ggplot()+geom_line( analysis_C,mapping=aes(x=scale(C_draw),y=yield),color='blue' )+geom_line( analysis_RUE,mapping=aes(x=scale(RUE_draw),y=yield),color='red' )+labs( title ="Sensitivity analysis",subtitle ="C (blue) and RUE (red)",x ="Parameters variation (relative scale)",y ="Potential max yield (t.ha-1)" )p3
Model optimization
Which parameter should be optimized first to improve model performance?